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Abstract. We present a study of the numerical solution of the two dimen- 
sional electrical impedance tomography problem, with noisy measurements of 
the Dirichlet to Neumann map. The inversion uses parametrizations of the 
conductivity on optimal grids. The grids are optimal in the sense that fi- 
nite volume discretizations on them give spectrally accurate approximations 
of the Dirichlet to Neumann map. The approximations are Dirichlet to Neu- 
mann maps of special resistor networks, that are uniquely recoverable from the 
measurements. Inversion on optimal grids has been proposed and analyzed re- 
cently, but the study of noise effects on the inversion has not been carried 
out. In this paper we present a numerical study of both the linearized and the 
nonlinear inverse problem. We take three different parametrizations of the un- 
known conductivity, with the same number of degrees of freedom. We obtain 
that the parametrization induced by the inversion on optimal grids is the most 
efficient of the three, because it gives the smallest standard deviation of the 
maximum a posteriori estimates of the conductivity, uniformly in the domain. 
For the nonlinear problem we compute the mean and variance of the maximum 
a posteriori estimates of the conductivity, on optimal grids. For small noise, 
we obtain that the estimates are unbiased and their variance is very close to 
the optimal one, given by the Cramer-Rao bound. For larger noise we use 
regularization and quantify the trade-off between reducing the variance and 
introducing bias in the solution. Both the full and partial measurement setups 
are considered. 



1. Introduction 

We study the inverse problem of electrical impedance tomography (EIT) in two 
dimensions, with noisy measurements of the Dirichlet to Neumann (DtN) map. 
Explicitly, we seek the positive and bounded, scalar valued coefficient a{x.) in the 
elliptic equation 

(1) V • [cr(x)V^/(x)] =0, X G 
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The domain Q is bounded, simply connected, with smooth boundary B. By the 
Riemann mapping theorem ah such domains in are conformahy equivalent, so 
from now on we take for Q the unit disk. We call cr(x) the conductivity and 
u e H^{Q) the potential, satisfying the boundary conditions 

(2) ii(x) = F(x), xgS, 

for arbitrary V G H^^'^{B). The data are finitely many noisy measurements of the 
DtN map : H^^^{B) H-^^^{B), which takes the boundary potential V to the 
normal boundary fiux (current) 

(3) A.y(x)=a(x)^|^, xeB. 

We consider both the full boundary setup, where A^rV^ is measured all around the 
boundary S, and the partial boundary setup, where the measurements are confined 
to an accessible subset Ba C S, and the remainder Bi = B\Ba of the boundary is 
assumed grounded {V\bj =0). 

In theory, full knowledge of the DtN map A^r determines uniquely a, as proved 
in ^34»d4| under some smoothness assumptions on a, and in [5] for bounded a. The 
result extends to the partial boundary setup, at least for a G e > as 

established in [24]. In practice, the difficulty lies in the exponential instability of 
EIT. It is shown in [Tl[3|33] that the best possible stability estimate is of logarithmic 
type. Thus even if the noisy data is consistent (i.e. in the set of DtN maps) we 
need exponentially small noise to get a conductivity that is close to the true one. 

It is shown in [2] that if a has finitely many degrees of freedom, more precisely 
if it is piecewise constant with a bounded number of unknown values, then the 
stability estimates on a are of Lipschitz type. However, it is not clear how the 
Lipschitz constant grows depending on the distribution of the unknowns in Q. For 
example, it should be much easier to determine the value of a near the boundary 
than in a small set in the interior of Q. 

An important question is how to find parametrizations of a that capture the 
trade-off between stability and resolution as we move away from the boundary, 
where the measurements are made. On one hand, the parametrizations should 
be sparse, with a small number of degrees of freedom. On the other hand, the 
parametrizations should be adaptively refined toward the boundary. 

Adaptive parametrizations for EIT have been proposed in [28l [31] and in [S] [4] . 
The first approaches use distinguishability grids that are defined with a linearization 
argument. The approach in [3l[4] is nonlinear and consists of an iterative coarsening 
and refinement of a piecewise constant discretization of the conductivity, with each 
discretization update being computationally costly. 

We follow the approach in [9l |36l [El [13] and parametrize a on optimal grids. 
The number of parameters is limited by the noise level in the measurements and 
their geometrical distribution in Q is determined as part of the inversion. The 
grids are based on rational approximations of the DtN map. We call them optimal 
because they give spectral accuracy of approximations of A^r with finite volume 
schemes. The grids turn out to be refined near the accessible boundary, where we 
make the measurements, and coarse away from it, thus capturing the expected loss 
of resolution of the reconstructions of a. 

Optimal grids were introduced in [6l [20l [211 [IS] for accurate approximations of 
the DtN map in forward problems. Inversion on optimal grids was first proposed 
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for Sturm-Liouville inverse spectral problems in [8 . The analysis in [11 shows that 
optimal grids are necessary and sufficient for convergence of solutions of discrete 
inverse spectral problems to the true solution of the continuum one. The numerical 
solution of EIT on optimal grids was introduced in P 123 for the full boundary 
measurements case, and in [121 [131 ES] for partial boundary measurements. The 
inversion in j9l [23l [121 [131 [32! based on the rigorous theory of discrete inverse 
problems for circular planar resistor networks [151 HSl [2S1 [HI US] , which gives net- 
works that can be uniquely determined by discrete measurements of the continuum 
DtN map [27l[9|. Just as in the continuum EIT, the inverse problem for networks is 
ill-posed, and there is a trade-off between the size of the network and the stability 
of the reconstruction. 

We present here a study of the inversion algorithms on optimal grids, for noisy 
measurements of the DtN map. We fix the number g of degrees of freedom, and 
analyze the effect of the adaptive parametrization of a on the reconstruction error. 
We consider maximum a posteriori estimates of a for the linearized problem about 
a constant conductivity, and for the nonlinear problem. The noise is mean zero 
Gaussian, and if its standard deviation is small, the only prior on a is that it is 
positive and bounded. For larger noise we use regularization (Gaussian priors), 
and study how the parametrization affects the trade-off between the stability of the 
result and the bias. 

We study three different parametrizations of a, with g degrees of freedom. The 
first two are piecewise linear, on an equidistant grid and on the optimal grid. The 
only relation between the second parametrization and resistor network inversion on 
optimal grids is the location of the grid nodes. The third parametrization is that 
induced by the resistor network inversion. 

In the linearization study we compute the standard deviation of the estimates 
and show that the resistor network parametrization is clearly superior. It gives 
estimates with uniformly small standard deviation in Q. The conclusion is that 
it is not enough to distribute the parameters on the optimal grid to obtain good 
results. To control the stability of the reconstructions, we also need to use proper 
basis functions. 

In the nonlinear statistical study we compute maximum a posteriori estimates 
of a with the inversion algorithms on optimal grids. We assess their quality by 
displaying pointwise in Q their mean and standard deviation. We obtain that 
the resistor network based inversion is efficient in the sense that it gives unbiased 
estimates of a, with variance that is very close to the optimal Cramer-Rao bound 
[35 . This is for small noise. For larger noise we use regularization priors that 
introduce some bias in the solution. We also compare the network based inversion 
to the usual optimization approach that seeks the conductivity as the least squares 
minimizer of the data misfit. For the optimization, the conductivity is piecewise 
linear with the same number g of degrees of freedom, either on a uniform grid or 
on the optimal grid. Our numerical experiments indicate that for a fixed allowed 
error (standard deviation) in the reconstructions, the network based method gives 
reconstructions that are closer in average to the true conductivity (i.e. with less 
bias). The conclusion for the non-linear problem is similar to that for the linearized 
problem: the reconstruction error is reduced with the network based inversion 
as compared to optimization on either equidistant or optimal grids. Our study 
considers both the full and partial measurement setups [9l [23l [121 [131 [32] • 
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The paper is organized as follows: We begin in section |2] with the estimation 
framework. Then we review the resistor network inversion in section |3l The tools 
needed for the numerical experiments are described in section [4j The estimation 
results are in section [5] We end with a summary in section |6] 

2. Maximum a posteriori estimation of the conductivity 

We study how different parametrizations of the unknown conductivity, with a 
fixed number g of degrees of freedom, affect the sensitivity of the reconstructions 
to noise in the data. Let s = (si, . . . , Sg)^ G be the vector of parameters, and 

(1) <7(X) = [5(s)] (x), XG^ 

the parametrization of the reconstruction conductivity, using an operator 5 : ^ 
C{Q) that takes s to a continuous function in Q. Since the data is noisy, the 
reconstructions are random variables. We study maximum a posteriori estimates of 
the reconstructions under certain priors, as explained in section (STTl We consider 



three different parametrizations of the form ([T]), outlined in section 2.2 



2.1. Estimation. We denote by F : C{Q) M^, the forward map that associates 
to a conductivity a the vector of measurements F(cr) of the DtN map. The measure- 
ment operation is explained in detail in section [5. 1.1 1 It amounts to recording volt- 
ages and currents at n electrodes on the boundary. The dimension g = n{n — l)/2 
of the data vector d corresponds to the number of independent measurements that 
can be made with n electrodes. The data model is 

(2) d = F(a) + e, eGAr(0,C), 

with € the noise vector. The notation A/'(0,C) states that e is Gaussian (multivari- 



ate normal), with mean zero and diagonal covariance C. We refer to section 3.1.1 
for an explanation of the uncorrelation of the components of e. 

We parametrize the conductivity as in ([T]) with a vector s G M^, with g = 
n{n — l)/2 being the same dimension as the data space. Since the data d is tainted 
with noise, we treat s as a continuum random variable, and denote by 7rp^(s) its 
prior probability density. The likelihood function 7r(d|s) is the probability density 
of d conditioned to knowing s. Given our Gaussian noise model, it takes the form 

(3) ^(d|s) = ^-^-l^exp [-i(F(5(s)) -d)^C-i(F(5(s)) -d) 

where \C\ is the determinant of the covariance C. The estimation of s is based on 
the conditional (posterior) density 7r(s|d). It is defined by Bayes' rule [551 

where 7r(s,d) is the joint probability density of (s,d). The marginal 

(5) 7r(d)= [ TT{s,d)TTjs)ds 

is just a normalization that plays no role in the estimation. The prior density tt^^ (s) 
may introduce a regularization in the inverse problem [29^ Chapter 3]. The priors 
used in our study are summarized in |A] They all ensure that the reconstructions 
(7(x) are positive. 
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We consider maximum a posteriori (MAP) estimates of the conductivity 

(6) Ct^apW = [5(S^AP)]W, 

which maximize the conditional probabihty density 7r(s|d). The vector s^^^p of 
parameters solves the optimization problem 

(7) s„,p = argmin[F(5(s)) - dJ^C"! [F(5(s)) - d] - log(7r^,(s)). 

The MAP estimates s^^^p are random, because they depend on the noise e in the 
measurements. To quantify their uncertainty, we approximate their variance using 
a large number M of independent samples s^^p , determined from ^ and data ^ 
with draws e^^) G A/'(0,C) of the noise, 

(8) Var[s] « ^ [s^ - (s>] , (s> « ^ E ^Lli- 

m=l m=l 

Then, we compare Var[s] to the optimal variance, which is the right hand side in 
the Cramer-Rao bound [35, Corollary 5.23], 

(9) E,.{[s„,p -E,.{s„,,}]2} > b^(s.)X-^(s.)b(s,). 

The notation Eg^ indicates that the mean (expectation) depends on the true vector 
of parameters s^. The bias factor b(s) is defined by 

(10) b(s) = VE3{s„,p}, 

where V denotes gradient with respect to s, and X(s) is the Fischer information 
matrix [33 Section 2.3]. It measures how much information the data d carry about 
the parameter s^. The Fisher matrix is in R^^^, with entries 

(11) X,,,(s) = Es {ds^ log^(d|s) ds^ log^(d|s)} . 

Since the likelihood 7r(d|s) is Gaussian, we obtain under the natural assumption 
that the noise covariance C is independent of s, that 

(12) X(s) = D,S{sYD^¥{S{s)fC-^D^¥{S{s))D,S{s), 

where DsS{s) is the Jacobian of S{s) evaluated at s and DfjY^cj) is the Jacobian 
of F(cr) evaluated at a. We take as the true the solution of the optimization 
problem ([7|), with noiseless data and upper /lower bound prior (22). The bias b(s^) 
is approximated via Monte Carlo simulations, from a large sample of draws. The 
bound ([9| is evaluated and displayed in section 5.2 



2.2. Parametrizations of the conductivity. We consider three different parametriza- 
tions: 

• piecewise linear on uniform grid: The entries in s are pointwise values 
of (j(x) on a uniform tensor product grid. The conductivity is piecewise 
linear on a Delaunay triangulation of the grid points. 

• piecewise linear on optimal grid: The conductivity and the parameters 
are defined as above, but with grid points from the so-called "optimal grid" 
(see section 3.2.1 ). 

• resistor network: The parameters s are (up to known multiplicative con- 
stants) the conductors in a network that has the same electrical response 
(DtN map) as the electrical response at the electrodes of the unknown con- 



ductivity. This parametrization is discussed in more detail in section 2.3 
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In the first two parametrizations the conductivity depends linearly on s. It is 
the piecewise linear interpolation of the entries in s, between grid nodes xi, . . . , x^, 

(13) ct(x) = [5(s)](x) = ^s,0,(x). 

i=l 

Here (pi are piecewise linear basis functions on a Delaunay triangulation of the nodes 
x^, satisfying the usual property 0i(xj) = dij^ with the Kronecker delta notation. 
In the resistor network parametrization the conductivity depends nonlinearly 



on the parameters s. The dependence is given in equation (16), in terms of the 
resistor network reduced model of ([T]). We emphasize that the resistor network 



parametrization cannot be written in the form (13). One way of comparing it to 
the linear case ( [l3| is to consider a small perturbation = {Ssi, . . . , Ssg)^ of some 
reference s, and linearize 

9 

(14) [S{s + Ss)] (x) = [S{s)] (x) + J2 ^^^^^ W + ^(^^)- 

Here the (pi are the "columns" of the Jacobian DsS{s) of the parametrization op- 
erator S with respect to s, and evaluated at s. We call these (pi sensitivity basis 
functions. 

2.3. The resistor network discretization and the sensitivity basis func- 
tions. For a given resistor network with n boundary nodes and g = n{n — l)/2 
resistors (one per edge), we denote by F : the discrete forward map 

which maps the vector of positive conductances of the network to a vector of g 
independent entries of the Dirichlet-to-Neumann map of the network. The choice 



and ordering of the independent entries is identical to that of F in section |2.1[ We 
assume that the topology (underlying graph) of the network is such that F admits 
a left inverse F~^, i.e. F~-'^(F(7)) = 7 for all 7 > 0. Such topologies are given in 
section I3.1.2[ 

In the resistor network parametrization we let 

s=(7i/7f\.--,7./7^'^ren. 

Here the 7^ are the conductances of the network that has the same electrical re- 
sponse as that measured for <j, i.e. F(<j) = F(7), for 7 = (71, . . . ,7^)"^. Similarly, 
7^"*-^ = (7^"''^ . . . , 7^"''^)-^ satisfies F(l) = F(7^-'-^). Hence, it is easy to compute s 
from knowledge of a and F~^. 

The mapping S in the resistor network parametrization is defined implicitly by 

(15) F(5(s))=F(s7«), 

where the product of the vectors s and understood componentwise. There 



are many functions that satisfy (15). We define S{s) for the resistor network 
parametrization as the limit of the Gauss-Newton sequence defined in (15), 



(16) [5(s)](x) = exp(/^(x)), where K, = \im K,j. 

The starting point in the Gauss-Newton iteration is k,o{^) = ln(<j^(x)), where 
cr^(x) is the piecewise linear interpolation of the discrete values s on the optimal 



grid defined in section 3.2.1 In practice, the evaluation of (|16[) is computationally 



efficient because iteration (15) converges quickly, basically after one step, as shown 
in [ills]. 
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Recall that for the resistor network parametrization, the mapping S is non- 



linear and thus cannot be written as an expansion in basis of functions, as in (13). 



To compare it to parametrizations of the form (13), we look at its linearization 



around reference parameters s to obtain (14). The basis functions (j)k are then the 



"columns" of DsS{s)^ which can be determined by differentiating (15) with respect 
to s, 

(17) D,S{s) = I),Ft(5(s))I)^F(s7^'^)diag(7^'^) = (diag(l/7^'^)I^a7)^ 



The sensitivity functions Dcr7(x) are defined in section 3.2.1 and are evaluated at 
a = S{s). They give the sensitivity of the resistors 7 to changes in the conductivity, 
and are an important ingredient of the inversion algorithm that is reviewed in 
section O 



2.4. Outline of the results. The linearized problem is studied in section 5.1 In 
this case it it is known [35, Section 5.1] that the Cramer- Rao bound is attained 
by the variance of the MAP estimates. The variance depends of course on the 
parametrization ([T]) used in the estimation. The results say that when we take 
piecewise linear interpolations of the parameters s on the equidistant or the optimal 
grids, we obtain much larger variances than if we use the sensitivity basis functions. 
That is to say, the sensitivity basis functions lead to better estimates than the 
piecewise linear ones, even when we interpolate on the optimal grids. We also show 



in section 5.2 that the nonlinear estimation based on resistor networks is efficient. 



in the sense that the sample variance is very close to the Cramer-Rao bound. 

3. EIT WITH RESISTOR NETWORKS 

The resistor networks used in our inversion algorithms are reduced models of ([T]) 
that are uniquely recoverable from discrete measurements of the continuum DtN 
map, as described in section |3.1[ They allow us to estimate the conductivity a 



parametrized on optimal grids, as explained in section 3.2 We give here a brief 



summary of the measurement operation, the critical resistor networks, and the 



induced reconstruction mapping ( 16 ) used in this paper. We refer to [9l [121 US] foi" 



more details of the resistor network based inversion. 

3.1. Resistor networks as reduced models for the forward and inverse 

problem. Resistor networks arise naturally in finite volume discretizations of equa- 
tion ([1]) on staggered grids with interlacing primary and dual grid lines that may be 
curvilinear. The potential is discretized at the primary nodes Pij^ the intersections 
of the primary grid lines, Uij ~ u{Pij). Each node Pij G Q is surrounded by a 
dual cell Cij, as shown in Figure [l] Integrating ([T| over the cells Q^-, using the 
divergence theorem, and approximating the boundary fluxes with finite differences 
we obtain a system of linear equations of the form 

Equations ([T]) are Kirchhoff's node law for the interior nodes of the resistor 
network with graph T = iV^S). Here V = {Pij} is the set of primary nodes, 
given by the union of the disjoint sets and V^- of boundary and interior nodes. 

Adjacent primary nodes are connected by edges, the elements of the set £ CV xV. 

I s I 

The network is the pair (r,7), with 7 G the vector with entries given by the 



(1) 
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Figure 1. Illustration of a staggered grid, with primary (solid) 
and dual (dashed) grid lines. The primary grid nodes are indicated 
by X and the dual grid nodes by o. We show a resistor as a 
rectangle with a midpoint □ at the intersection of a primary and 
dual line. 



conductances (inverse of resistances) ja,i3 > of the edges, following a preassigned 
ordering of S. Here (a,/3) ^ {(^, j ± ^) , {^^IJ)}' 
We relate the conductances ^^,(3 and a by 

(2) 7«,/3 = <T(P„,,)^^, 

where L denotes the arclengths of the primary edges Ea^i3 and dual edges T^a,/3' The 
points Pa, 13 are located at the intersections of the primary and dual grid segments 

Ea,(3 and ^a,(3' 

The forward problem for a known network (r,7) amounts to determining the 
potential function U : V ^ with Uij = V({Pij) satisfying the conservation of 
currents ([T]) at the interior nodes, and Dirichlet boundary conditions 

(3) u\rs=UB' 

We denote the number of boundary nodes by n. The entries in the vector G 
may be related to the continuum boundary potential V as explained below, in 
section I3.1.1[ 

The inverse problem for the network seeks the conductances 7 from the discrete 
DtN map A^. The graph T is known, and the DtN map is a matrix in R^^^ 
that maps the vector Uj^ of boundary potentials to the vector of boundary 
current fluxes. Since we consider the two dimensional problem, all the graphs T are 
circular planar graphs [15l[T6], i.e. graphs that can be embedded in the plane with 
no crossing edges and with all boundary nodes Vb lying on a circle. 

3.1.1. Discrete measurements of the continuum DtN map. To connect the 
discrete inverse problem for the network (r,7) to continuum EIT, we introduce 
a measurement operator Mn that defines a matrix Mn (Act) G W^^^ from the 
continuum DtN map A^r. The measurement operator is chosen so that for any 
suitable conductivity, M.^ (^cr) is consistent with the DtN map of a circular planar 
resistor network (r,7). The network is a reduced model for the forward problem. 
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because it satisfies 

(4) = Mn (A.) . 

The continuum forward map F of section [2] is defined using this measurement 
operator as 

(5) F{cT)=vec{Mn{Aa)). 

where Yec{A) denotes the operation of stacking in a vector in R^, ^ = n(n — l)/2, 
the entries in the strict upper triangular part of a matrix A G R"^^"^. Because 
of reciprocity and conservation of currents, these entries completely determine the 
measured DtN map A4n (Act). Hence, another (equivalent) way of writing the 
compatibility condition Q is 

(6) F(7) = F(a). 

where F(7) = vec(A^) is the discrete forward map. 

One possible choice of the measurement operator consists of taking point values 
of the kernel of A^r. Its consistency with networks is shown in [27\ [25]. Another 
choice, which we use in this paper, is to lump fluxes over disjoint segments of B 
that model electrode supports. Its consistency with networks is shown in [9l [23]. 
Such an operator is defined using the nonnegative "electrode" functions xi, • • Xn 
in with disjoint supports, numbered in circular order on B. We normalize 

them to integrate to one on B. The operator A4ri maps A^r to the symmetric matrix 
with off-diagonal entries given by 

where (•, •) is the duality pairing between H^^'^{B) and R-^^'^^B). The diagonal 
entries are taken so that the rows (and columns) of A^n(Acr) sum to zero. Such 
choice enforces the conservation of currents. 

3.1.2. Solvability of the inverse problem for resistor networks. The ques- 
tion of solvability of the inverse problem for circular planar networks like (r,7) has 
been settled in [151 [SI [23 [H [H] • The answer is that when the graph T is critical^ 
the discrete forward map IF is one-to-one, and there exists a left inverse so that 
F~-^(F(7)) = 7 for all 7 > 0. A graph is critical if it is well connected and if it does 
not contain any redundant edges. See jT5l[T6] for a technical definition of criticality 
and well-connectedness. In a critical network, the number n of boundary nodes and 
the number g of edges in the graph obey g = n{n — l)/2. This says that there are 
as many unknown conductances in the network as there are degrees of freedom in 
the DtN map A^ G R^><^. 

It remains to define the graph of the network, so that it is critical, and thus 
uniquely determined by Q. Typically, different graph topologies are better suited 
for the full and partial data measurements. Here we use the topologies considered 
in ^ |23l [la [32l |10 , see [B] for details. 

The inverse problem for critical networks (r,7) can be solved with at least two 
approaches. We use them both in our study of the nonlinear inverse problem in 
section |5l 

(1) Layer peeling [151 |13l HO] • A direct method giving the conductances 7 
in a finite number of algebraic operations. The advantage of layer peeling 
is that it is fast and explicit. The disadvantage is that it becomes quickly 
unstable, as the size of the network grows. Moreover, noisy data may not 
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be consistent with a network, i.e. the consistency relation ([6| may not hold 
if the forward map F is noisy. 
(2) Optimization: Use standard optimization techniques to find conductances 
7 that best fit the (possibly noisy) measurements in the least squares sense 



(see section 4.2 for more details). 



3.2. Inversion on optimal grids. We denote by the true conductivity, to dis- 
tinguish it from the estimates that we denote generically by a. The relations ([2| 
between 7 and have been derived in the discretization of the forward problem. 
We use them for the conductances of the network (r,7) recovered from the mea- 
surements F(cr^), in order to estimate a^. This does not work unless we use a 
special grid in ([2| [9^ 11 . The idea behind the inversion on optimal grids is that 
the geometrical factors L(Sq,^^)/L(£^q,^^) and the distribution of the points P^^is in 
([2]) depend weakly on a. Therefore, we can determine both the geometrical factors 
and the grid nodes from the resistor network (r,7'^^^), with the same graph F as 
before, and F(7(i)) = F(l). These are the measurements of the DtN map for con- 
stant conductivity a = 1, that we can compute, and 7^^^ ^ L{T^a,/3)/L{Ec^i3). We 
obtain the pointwise estimates 

(8) ^(Pc.,/3) 



that we place in Q at points Pck,/3 determined from a sensitivity analysis of the DtN 
map, as we explain next. 

3.2.1. The sensitivity functions and the optimal grids. The distribution of 
points Pa, (3 in Q is optimal in the sense that 

(9) F(7(ct)) = F(ct) 

for conductances 7 = 7(cr) related to the continuum a as in ([2]), and for a = 1 (i.e., 
for 7(1) = 7*^^^). Each conductance is associated with a point Pa,/3^ so we write 
k = k{a^/3). We define the optimal grid points as the maxima of the sensitivity 
functions given below, evaluated at a = 1, 



(10) ^«,/3 = argmax(L>^7^(^ ^))(x) 



El 



These are the points at which the conductances are most sensitive to changes in 
the conductivity. 

To compute the sensitivity functions, we take derivatives in ([9| with respect to 
cr, and obtain 

(11) = {D^W{f{a)))-' D,F{a). 

The left hand side is a vector-function from (7 to R^. Its k—th entry is the sensitivity 
of conductance 7^ with respect to changes of a. The matrix D^¥{'y{(j)) G M^^^ 
is invertible [15]. The Jacobian DcrF{a) can be written in terms of the Green's 
function of the differential operator u V-(crVii) with u\jg =0, and the "electrode" 



functions Xi introduced in section 3.1.1 The calculation is given in detail in J3l 
Section 4]. 
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3.2.2. The estimate of the conductivity on optimal grids. What we have 
computed so far allows us to define an initial estimate cr^(x) of the conductivity, 
as the linear interpolation of the values ([8|, on the optimal grid defined by (10). 



Then, we improve the estimate using a Gauss-Newton iteration that minimizes the 
objective function 

(12) J{a) = mF{a))-Q{d)\\l 

over search conductivity functions cr(x). Here the data vector is given as in (H for 
the conductivity cr^, that we wish to find. The reconstruction mapping is 

(13) Q(d) :=diag(l/7(i))F-i(d), 

and it involves solving the discrete inverse problem for a resistor network. The 
map Q computes the pointwise estimates (|8| from the data d. Therefore, cr^(x) is 
a linear interpolation of Q(d) on the optimal grid. 

The objective function J{cr) is different than the usual output least squares data 



misfit ||F(cr) -d||2. We use Q in p2l as a nonlinear preconditioner of the forward 



map F, as explained in detail in |9l|23J. It is because of this preconditioning, and 
the good initial guess cr^(x), that we can obtain close estimates of a by minimizing 
J7(cr) in [9l[23] and in this paper. The estimates are computed with a Gauss-Newton 
iteration that basically converges in one step [Qj [23] . 

We enforce the positivity of a by the change of variables = ln((j), so that we 
work with the map 

(14) G(/t) = ln[Q(F(exp(K)))]. 



The Gauss- Newton iteration that we use in the parametrization (16) is 
(15) = Kj^i + {D^G{kj^i))^ [ln[Q(d)] - G{kj^i)] , j = 1, 2, . . . , 

with initial guess = In(cr^). The index f denotes the Moore-Penrose pseudo- 
inverse, and the iteration amounts to finding the update hvj K/j — \ as the orthogonal 
projection of the residual onto the span of the sensitivities, the column space of the 
transpose of D^G{hij-i). These sensitivities are easily related to those computed 



in section 3.2.1, using the chain rule to deal with the change of variables = In(cr). 
4. Numerical experiments setup 



We explain in section [471] how we simulate the noisy measurements. The noise 
may be too high for the layer peeling method to work. The optimization method 
presented in section [4.2[ is more robust to noise and, as a bonus, it allows us to 
solve efficiently the optimization problem for the MAP estimate ^ with the re- 
sistor network discretization (see Remark [T]). For reference, we include noiseless 
reconstructions in section [431 



4.1. Data and noise models. We solve equation ([T]), with a second order finite 
volume method on a very fine, uniform, tensor product grid, with N n nodes 
on the boundary. This approximates A^Ar(Acr^), where A^at is the measurement 
operator of section 3.1.1 The noise e in ([2| is given by 

(16) e = vec {Mn{S)) , S e R^""^. 

We use two noise models, defined in terms of the noise level £ and a symmetric 
N X N matrix 77, with Gaussian, identically distributed entries with mean zero and 
variance one. The entries in 77 on and above the diagonal are uncorrelated. The 
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first model scales the noise by the entries of the DtN map for constant conductivity 

(7 = 1, 

(17) S=iMN{Ai)-V. 

where symbol • stands for componentwise multiplication. The scaling makes the 
noise easier to deal with, and the model is somewhat similar to multiplicative noise. 
The second noise model is 

where || || is a matrix norm that approximates the continuum H^^'^(B) H~^^^{B) 
operator norm. It is defined in[Cl 

4.2. The Gauss-Newton iteration for determining the resistor networks. 

The direct, layer peeling algorithms described in [151 [131 ttO] are fast, but highly 
unstable and can be used only for very small noise. It is not known how to regularize 
layer peeling algorithms. To deal with the instability, we can only reduce the size of 
the network, as was done in [9]. However, simply reducing the network size is not 
sufficient for the larger noise levels considered in the simulations. We use instead the 
more robust Gauss-Newton method described below, which allows regularization. 
The Gauss-Newton method is more expensive than layer peeling (about 20 times 
more expensive for n = 7), but the computational cost is reasonable because the 
dimension g of the vector of unknown conductances is small for noisy data, and the 
Jacobian is relatively inexpensive to compute. 

The Gauss-Newton method determines the log-conductances n = ln(7) in the 
network (r,7) (with topology F fixed) by minimizing the objective functional 

(19) 0{k) = (F(exp(A^)) - d)^C-^(F(exp(A^)) - d) + a| — A^ref 1 1 7 

over K, e R^. Thus, the positivity of the conductances 7 = exp(A^) is satisfied 



automatically. The first term in (19) measures the misfit between the measured 
data d modeled by ([2|, and the data produced by a network with log-conductances 



K. Here C is the covariance matrix of the measurements. The second term in (19) 
is a Tikhonov type regularization penalizing the distance from some reference log- 
conductances A^ref- The parameter a > determines the strength of the penalty 
term. 

In our numerical experiments, the Gauss-Newton approximation of the Hessian 
of 0{k) is further regularized by adding 10~^/Ci^i to its diagonal. The iterations 
are stopped either when the norm of the gradient of 0{k) is 10~^ smaller than that 
at the initial iterate, or when the maximum number of iterations (300) is reached. 
The initial iterate is k, = ln(7^^^), and we take as reference log-conductances K^-ei = 
ln(7(i)). 



Remark 1. Recall from (15) that the network is the reduced model that matches 
the data. Therefore: 

(1) When a > 0; finding log- conductances that minimize 0{k) is equivalent to 
finding the MAP estimate with the resistor network based parametriza- 



tion, and the prior ( f^^p . The regularization parameter a is the same in 
both ^ and (19). Moreover, when Kref = In (7^-^^) in (19), the reference 
parameters appearing in (24) are Sref = 1 (vector of all ones with length g). 
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0.5 1 1.5 2 0.5 1 1.5 2 2.5 

Figure 2. Conductivities used in the numerical experiments. 
Left: smooth conductivity. Right: piecewise constant chest phan- 
tom. 



(2) When a = 0, minimizing 0{k,) is equivalent to finding the MAP estimate 
^ with the resistor network based parametrization and the upper/lower 
bound (for a large enough upper bound ama^)- 

4.3. Noiseless reconstructions. We show in Figures [3] and [4] the reconstructions 
of the true smooth and piecewise constant conductivities displayed in Figure |2j 
The reconstructions are obtained in the full and partial boundary setups, and for 
noiseless data. Two distinct cases of partial boundary measurements are considered. 
In the one-sided case the accessible boundary Ba consists of a single connected 
segment of B. In the two-sided case the accessible boundary Ba consists of two 
disjoint segments of B. 

We display in the top row of Figures [s] and [4] the initial guess cr^(x) of the 



Gauss-Newton iteration (15). It is the piecewise linear interpolation of the values 



o'^{Pa,p) = <^k{a,i3)^ where cT/c are obtained from and the optimal grid nodes 



Pa^^ are defined by (10). The function cr^(x) is linear on the triangles obtained by 
a Delaunay triangulation of the points Pa, (3- In the partial boundary measurements 
case we display cr^(x) in the subdomain delimited by the accessible boundary and 
the segmented arc connecting the innermost grid points. We set cr^(x) to the 
constant value one in the remainder of the domain. The plots in the bottom row 
in Figures [3] and |4] display the result of one step of the Gauss-Newton iteration 
described in section [3.2.21 

The network topologies are defined in B] The reconstructions from full boundary 
data are in the left column in Figures |3 and[4j They are obtained with a circular 
network C (^^^,n), with n = 29. The reconstructions with the one-sided partial 
boundary measurements are in the middle column, and they are obtained with a 
pyramidal resistor network F^, for n = 16. The reconstructions with the two-sided 
boundary measurements are in the right column. They are obtained with a resistor 
network T^, for n = 16. 

5. Effect of parametrization on the reconstruction error 

We study numerically the effect of noise on the inversion with (a) resistor net- 
works, and (b) the conductivity parametrized with piecewise linear basis functions. 
We consider in section [STT] the linearized problem about the constant conductivity 
a = 1, and in section [5^2] the nonlinear problem. 
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Figure 3. Reconstructions of the smooth conductivity with 
noiseless data. 




Figure 4. Reconstructions of the piecewise constant conduc- 
tivity with noiseless data. Here and in Figure [3j Left column: full 
boundary measurements. Middle column: one-sided partial mea- 
surements. Right column: two-sided partial boundary measure- 
ments. The top row shows the initial guess cr^(x). The bottom 
row shows the result of one step of the Gauss-Newton iteration. 
The color scale is the same used for the true conductivity in Figure 



5.1. The linearized problem. The results in this section are for linearization at 



the constant conductivity a(x) = 1, with additive noise modeled as in (18). All 
the three parametrizations described in section 2.2 represent exactly a = 1, with 
a parameter vector s of all ones, i.e. S{1) = 1. Hence, the linearization of the 
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Figure 5. Standard deviation of reconstructions for the lin- 
earized problem, (a): Piecewise linear basis functions on uniform 
grid, (b): Piecewise linear basis functions on optimal grid, (c): 
Resistor network parametrization. The color scale is base-10 log- 
arithmic and the noise level is ^ = 0.01% (additive noise model 
(18)). The circular resistor network has n = 11 boundary points 
and the noise E is given at A/" = 100 equidistant points on the 
boundary. 



forward map around s = 1 can be written as 

(20) F(5(l + (5s)) = F(l) + I),F(l)I)s5(l)(5s + o((5s), 

for a small perturbation of the parameters. If we discretize the conductivity on a 
fine grid with A^^r points, the Jacobian D^Y{V) is a x A^^r matrix, and the Jacobian 
D^S(X) is an A^^r x Q matrix with columns given by the basis functions in ([l4|. 
The estimate cr^^p is calculated by solving the optimization problem ([7|), with 



the linearization (20) of the forward map, and the upper /lower bound prior ( |22| ) 
This optimization is a quadratic programming problem that we solve using the 
software QPC [37]. The mean and variance of cTj^ap estimated with Monte 
Carlo simulations 

(2¥).r[a(x)] « ^ [a^ (x) - (<x(x))] , (a(x)) « - ^ (x), 

m=l m—X 

using M = 1000 samples. We do not show the mean (cr(x)) because it is basically 

1/2 

a(x) for all the cases that we present below. The standard deviation {Var [cr(x)]} ^ 
is shown in Figure [5] for the case of full boundary measurements, and the three 



parametrizations described in section 2.2 The noise level is ^ = 0.01%. We choose 



it so small to minimize the action of the positivity constraints imposed by the prior. 
Larger noise levels are considered later in the paper. 

With the resistor network parametrization, the standard deviation is smaller and 
does not increase toward the center of the domain. Also there are no active positiv- 
ity constraints. The random fluctuations of cTj^ap Figure [5](b) lie mostly within 
three standard deviations, and are much smaller than the background conductivity 
a = 1. 

The piecewise linear parametrization on the equidistant grid gives a large, order 
one standard deviation in the center of the domain. The positivity constraints are 
active in 81.6% of realizations. Surprisingly, the piecewise linear parametrization 
on the optimal grid is worse. Its standard deviation is large, of order one in most 
of the domain, and the positivity constraints are active in 61.1% realizations. This 
shows that it is not enough to distribute the parameters on the optimal grid. 
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Figure 6. Comparison of the condition number of the matrix 
DcrF{l)DsS{l). The abscissa is the number n of boundary nodes 
and the ordinate is in logarithmic scale. The green line is for 
the fine grid forward map {S = identity). The red line is for the 
piecewise linear functions on the equidistant grid. The blue line is 
for the piecewise linear functions on the optimal grid. The black 
line is for the resistor network parametrization. 




Figure 7. The standard deviation (base- 10 log scale) for one- 
sided (left) and two-sided (right) boundary measurements. The 
resistor networks have n = 16 boundary points and the noise S 
with level £ = 10'^% is given at TV = 83 (left) and N = 104 
(right) equidistant points on the accessible boundary. Top: piece- 
wise linear parametrization on the optimal grid. Bottom: resistor 
network parametrization. The accessible boundary is in solid red. 

The same conclusion can be reached from Figure [6} where we display the con- 
dition number of the matrix DcrF{l)DsS{l), as a function of the number n of 
boundary points. The condition number increases exponentially with n, as ex- 
pected from the exponential ill-posedness of the problem. However, the rate of 
increase is smaller for the resistor network parametrization. 

The standard deviation {Var[cr(x)]}^^^ is shown in Figure [t] for the case of one 
and two-sided boundary measurements. We use a much smaller noise level {i = 
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10~ %, additive model (18)) in the partial measurements case than in the full 
data case, because we compute the standard deviation for bigger networks (n = 
16, same network size as in the reconstructions of Figures |3]and[4|. Since the 
condition number of the linearized problem grows exponentially as we increase n, 
only very small levels of noise can be used for n = 16. We present the results for 
the piecewise linear parametrization on the optimal grid (top row) and the resistor 
network parametrization (bottom row). 

We reach the same conclusion as before. The resistor network parametrization 
gives a smaller standard deviation, that does not increase toward the inaccessible 
region. The piecewise linear parametrization on the optimal grid gives a large 
standard deviation near the inaccessible region in the one-sided case, and in the 
whole domain in the two-sided case. The positivity constraints are active in most 
realizations for the piecewise linear parametrization. They are not active for the 
resistor network parametrization. 

5.2. The nonlinear problem. We study the statistics (mean and standard de- 
viation) of the MAP estimates of the parameters s^^p and the conductivities 
^MAP = ^^(Smap)^ which come from minimizing the functional ([7|. The study can 
be done for both the full and partial boundary measurement setup, but we present 
here only the full measurements case. We consider first, in section [5.2. 1[ very small 
noise so that we can use the fast layer peeling inversion algorithm to find the min- 
imizer of ([7|), with resistor network parametrization and upper/lower bound prior 



(22 ). The speed of the algorithm allows us to compute the Cramer- Rao lower bound 
in a reasonable amount of time. We do not calculate the bound for larger noise, 
where we use regularized Gauss- Newton to determine the resistors, because of the 



computational cost. However, we do show in section 5.2.2 the bias and relative 



standard deviation of the reconstructions, and we also compare in section 5.2.3 the 



results to those with a piecewise linear discretization and Gaussian prior on the 



conductivity (23). 



5.2.1. Statistics of resistor network inversion using layer peeling. The 

results in this section are for the MAP estimates s^^p, solving the optimization 
problem ([t]), with the resistor network parametrization, and upper/lower bound 
prior ( [22| ). We present the mean (s), Bias[s] = — (s), and the standard deviation 
(Var[s])^/^ of the estimates. We consider a very small noise level £ = 0.1% (noise 



model (p^), so that we are able to minimize (19) with no regularization {a = 0) 
directly, using the layer peeling algorithm. 



Note that in the Cramer-Rao bound, the Fischer matrix (11) can be calculated 



analytically, but the bias factor (10) is estimated with Monte Carlo simulations. 
This is the expensive part of the computation, because we need a large number 
of samples to estimate the mean. In addition, each component of the vector s is 



perturbed to approximate the partial derivatives in (10) via finite differences. We 
use M = 1000 samples, and the bias is relative to s^, the solution of the optimization 
problem ^ with noiseless data. The partial derivatives in (10) are approximated 
with finite differences with a step size of 0.01. 

The bias factor is shown in Figure |8| and it is close to the identity matrix. That is 
to say, the estimates are unbiased. Figure^ shows (a) the mean (s), (b) Bias[s] and 
(c) the relative standard deviation (Var[s])^/^/ (s), where the division is understood 
componentwise. The last column (d) shows the difference in percentage between 



18 



L. BORCEA, F. GUEVARA VASQUEZ AND A.V. MAMONOV 




piecewise constant 



Figure 8. Estimated bias factor (10) for conductivities in 
Figure [2] is close to the identity. 





lis 



2 3 10-4 lo'^r^3io'-2io'-^oo 2% 4^^o 

(b) (c) (d) 



Figure 9. Results with the direct resistor finding algorithm: 
(a) mean, (b) absolute bias, (c) standard deviation relative to the 
mean, (d) relative difference between variance and Cramer-Rao 
bound relative to Cramer- Rao bound (in percentage). The noise 
model is (17) and the level is ^ = 0.1%. The circular resistor 



network has n = 11 boundary points and the noise matrix S is 
given at = 100 equidistant points on B. All parameters are 
linearly interpolated on the optimal grid. 



Var[s] and the Cramer- Rao bound in ([9|, normalized pointwise by the Cramer- Rao 
bound. We evaluate the Cramer- Rao bound by setting the bias factor (10) to the 
identity, which is a good approximation (recall Figure [8| . Note that the difference 
between the variance and the Cramer-Rao bound is very small, indicating that the 
estimation is efficient. The result in column (d) should be non-negative. We have 
some negative numbers, probably due to insufficient sampling in Monte Carlo, but 
they are so small in absolute value that we can treat them as essentially zero. 



5.2.2. Statistics of conductivity estimates using optimization. Here we con- 
sider the additive noise model (18), with noise levels i = 0.01% and i = 1%. These 
are the levels used in [30 , and we use them to compare our results with those in 
[30]. Because solving ^ with only the upper/lower bound prior (22) does not give 
reliable estimates, we also use the prior (24). By Remark [l] this is equivalent to 
minimizing (19), which is computationally cheaper. 
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0.01% additive noise (18) and a 
^ in (19). The circular resistor 



Figure 10. Results for £ 
regularization parameter a = 10 
network has n = 11 boundary points and S is given at = 100 
equidistant points on B. Top: the smooth conductivity shown on 
the left in Figure [2] Bottom: the piecewise constant conductivity 
shown on the right in Figure |2j (a) mean, (b) bias, (c) relative 
standard deviation. 
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: 1% additive noise (18) and a 
in (19). The circular resistor 



Figure 11. Results for £ = 
regularization parameter a = 10 
network has n = 7 boundary points and 8 is given at = 100 
equidistant points on B. Top: the smooth conductivity shown on 
the left in Figure [2] Bottom: the piecewise constant conductivity 
shown on the right in Figure [2] (a) mean, (b) bias, (c) relative 
standard deviation. 



Figure 10 shows (a) the mean ((t(x)), (b) Bias [a (x)] and (c) the relative stan- 
dard deviation (Var[cr(x)])^/^/ (cr(x)), for the noise level £ = 0.01%. The bias is 
computed with respect to a^, the solution of the optimization problem ([7|, with 
noiseless data, no regularization {a = 0), and a resistor network parametrization 
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Figure 12. Realizations of cTj^apC^) for n = 7 boundary points, 
i = 1% additive noise (18) and regularization parameter a = 10~^ 



in (19). 



with n = 11 boundary nodes. The regularization parameter a = 10~^ is chosen 
so that both the bias and the relative standard deviation are small. The choice of 
the regularization parameter is discussed in more detail in section [5. 2. 3[ We do not 
show realizations of the MAP estimates, because they are close to the mean, as the 
standard deviation of the reconstructions is below 10~^. 



For the higher noise level i = 1%, we present in Figure 11 the results for a 
smaller network, with n = 7 boundary nodes. Again we choose the regularization 
parameter a = 10~^ in such a way that both the bias and standard deviation are 
small. The relative standard deviation is less than 10%, and the realizations of 
^MAp(^) shown in Figure 12 resemble the mean in Figure 11 These realizations 
are comparable to the reconstructions in [30 . The reconstructions with n = 11, 
noise level £ = 1% and an appropriate choice of the regularization parameter are 
qualitatively similar to those with n = 7 shown in Figure [TT] and thus are not 
included here. 

5.2.3. Resistor network parametrization compared to other parametriza- 
tions. We now study the interplay between regularization and parametrization. 
We solve the optimization problem ^ with the three parametrizations of sec- 
tion 12.21 At the noise levels considered here, the same as in section 15.2.2 



we 



need regularization to get reliable estimates of the conductivity. As was the case 
in section [5. 2. 2[ the reconstructions with the resistor network parametrization are 
regularized with the prior (24). For the piecewise linear parametrizations, we regu- 
larize with the Gaussian prior (23), with reference conductivity <Jref = 1. Moreover, 



we stopped the iterations when either the maximum number of iterations (70) is 
reached, or when the norm of the gradient of the objective function at the current 
iterate is smaller by a factor of 10~^ than that at the initial iterate. We also added 
10~^/Ci,i to the diagonal of the Gauss-Newton approximation to the Hessian. Re- 
call that C is the covariance of the measurements. 

We use two metrics to evaluate the reconstructions using different parametriza- 
tions. The first one is the norm of the true bias 

(a(x)))2dx''^' 



TrueBias[cr] = (^J (cr^(x) 



as a percent of the true conductivity cr^. This measures the fidelity (in average) of 
our reconstructions. The second metric is the norm of the standard deviation 
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Figure 13. The true bias (in abscissa and in percent) against 
the relative standard deviation (in ordinate and in percent) for 
different values of the regularization parameter. The curves corre- 
spond to the resistor network approach (black), the linear interpo- 
lation on the optimal grid (blue) and the uniform grid (red). The 
number of realizations N = 100. The regularization parameter 
was a = 10"-^/^, with j = 1,2, ... ,6 for the linear interpolation 
(red and blue) and with j = 1, 2, . . . , 12 for the network approach 
(black) . The regularization parameter used for the reconstructions 



in figures 10 and 11 is indicated with a black star. 



relative to norm of the mean of the reconstructions 

^ ^ (LVar[a(x)]dx)'/' 
RelStd[(j] - L V ;j y 



(L(^(x)) 



1/2 



which measures the stability of our reconstructions. 

We report in Figure [T3| for different values of the regularization parameter a, 
the true bias versus the relative standard deviation. Note the typical L-curve 
shape, which refiects the trade-off between accuracy (small bias) and stability (small 
standard deviation). When the regularization is not sufficient, the bias is small 
but the standard deviation is large (vertical branch). When the problem is over 
regularized, the bias is large and the standard deviation is small (horizontal branch). 
The "best" choice of the regularization parameter would be near the "corner" of 
the L shape. 

The first row in Figure 13 shows that for noise level 1 = 1% and for n = 7 
boundary nodes, the resistor network parametrization outperforms the piecewise 
linear parametrizations: for a fixed standard deviation, the bias is smaller. This is 
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specially noticeable for the piecewise constant conductivity. Interestingly, the uni- 
form grid is slightly better than the optimal grid in this case. For n = 11 boundary 
nodes and I = 0.01% (second row), all approaches give comparable results, with 
the resistor network giving bias smaller by a few percent, specially if we allow a 
standard deviation above 0.1%. 

Remark 2. We emphasize that the cost of solving ^ with the resistor network ap- 
proach (i.e. solving JT^^ is negligible compared to computing the Jacobian DsS{s). 
Thus, the resistor network approach takes about the same time as one step of Gauss- 
Newton to solve ^ with the piecewise linear parametrizations. In the computations 
for Figure [7^ the mean number of iterations for these linear parametrizations was 
at least 8, and varied depending on the regularization parameter, the grid and the 
conductivity. Therefore, the resistor network method is at least 8 times faster than 
the one using conventional discretization. 

6. Summary 

We presented a numerical study of the effects of noise on resistor based inversion 
algorithms. The algorithms were introduced in [H |23l [TJl [131 132] , and are briefly re- 
viewed here. We have three measurements setups. The first assumes that the entire 
boundary B of the domain Vt is accessible. The other two are for partial boundary 
measurements confined to the accessible boundary Ba C B. One setup assumes 
one sided measurements, with Ba consisting of a segment of B. The inversion al- 
gorithm is introduced in [131 IS2] • The other setup is two sided, with Ba consisting 
of two disjoint segments of B. The inversion amounts to defining a reconstruction 
mapping 5(s), that takes a vector s = 7/7^^^ of ratios of positive conductances 
of a network (r,7) and a reference network (r,7*^^^), to continuous conductivity 
functions defined in Vt. The network has a special graph V that is adapted to the 
measurement setup and which allows the conductors 7 to be determined uniquely 
from measurements of the DtN map. The mapping 5(s) involves a Gauss- Newton 
iteration that minimizes a preconditioned data misfit in the least-squares sense. 

Our study considers three different parametrizations of the unknown conductiv- 
ity with g degrees of freedom. The first two are piecewise linear interpolations on an 
equidistant grid and on the optimal grid, respectively. The third parametrization 
is based on resistor networks. 

For the linearized problem, the piecewise linear parametrizations give large vari- 
ances of the MAP estimates, even if we use the optimal grids. The resistor network 
parametrization is superior because the variances of the MAP estimates are lower 
and do not increase toward the inaccessible part of the domain. 

The statistical study of the non-linear problem shows that when no additional 
prior (regularization) is introduced, and the noise is very small, using the resistor 
network parametrization gives reconstructions with small bias and the variance of 
the MAP estimates is very close to the optimal Cramer-Rao bound. For larger 
noise, we regularize the resistor based inversion with a prior on the conductances, 
and we compare the results with those of output least squares with piecewise linear 
parametrizations of the conductivity, on uniform and optimal grids and regularized 
with a Gaussian prior on the conductivity. The study assumes realistic noise levels 
[30]. All three parametrizations give a trade-off between accuracy (small bias) 
and statistical stability (small standard deviation of the estimates). However, the 
resistor based parametrization consistently outperforms the piecewise linear ones. 
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giving a smaller bias for a fixed standard deviation in the reconstructions. The 
quality of the reconstructions is comparable to that in [30 . 

Our regularization priors are very simple. If additional prior information is 
available, the results of the resistor based inversion can be greatly improved, by 
enlarging the space of the Gauss-Newton iterates, beyond the span of the sensitivity 
functions of the reconstruction mapping. This was shown e.g. in [9, §7.1]. 

From the computational point of view, the inversion with resistor networks can 
be done at roughly the cost of one Gauss-Newton iteration for a conventional out- 
put least squares method, with the same number of degrees of freedom of the 
parametrization. In our numerical experiments, the resistor network inversion was 
at least eight times faster. 
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Appendix A. Prior distributions 

(1) Upper/lower bound prior: We use this prior alone to explore the effect 
of the parametrization on the stability of the reconstructions, at small noise 
levels. It states that the conductivity is positive and bounded. Let S = 
{s G : [5(s)](x) G (0, (Jmax), X G ^1} be the set of parameters mapped 
by S to positive conductivity functions bounded by (Jmax in ^- The prior 
is 

(22) "rHs) = ^, 

where ls(s) is the indicator function that takes value one when s G S, 
and zero otherwise, and |S| is the volume of S. When we study maximum a 
posterior estimates of the conductivity in section [2]T] we set (Jmax to a large 
enough value, and keep at the same time the number g of parameters low 
enough, for the constraint [<S(s)](x) < (Jmax to be automatically satisfied. 
However, we do enforce the posit ivity. 

(2) Gaussian prior on the conductivity: This is a Tikhonov regularization 
prior that is useful at higher noise levels [29, Chapter 3]. It says that in 
addition to the conductivity being positive and bounded, we assume that a 
has a normal distribution with mean (Tref. The fiuctuations S{s) — (Jref are 
uncorrelated from point to point, and the pointwise variance is a~^. The 
prior is defined by 

(23) 7r(f^)(s) ~ Trf ^)(s)exp[-(a/2)||5(s) - <T,ef 

where the symbol means equality up to a positive, multiplicative con- 
stant. 
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Figure 14. Circular resistor networks C(/,n) with critical 
graphs: I = {n — l)/2. The interior nodes are indicated with 
dots and the boundary nodes with crosses. 



(3) Prior on the parameters: This is also a Tikhonov type regularization 
prior that is useful at higher noise levels. It says that in addition to the 
conductivity being positive and bounded, the vector of the logarithm of the 
parameters s is normally distributed, with mean log(sref) and covariance 
a~^I, where / is the g x g identity matrix, 

(24) 7r(f^)(s) ~ Trf ^)(s)exp[-(a/2)|| log(s) - log(s,ef)||l]. 



Appendix B. Resistor network topologies 

Resistor networks (r,7) with circular graphs T = C{l^n) are natural reduced 
models of the problem with full boundary measurements. The notation Cil^n) 
[ini [17] indicates that the graph has / layers, and n edges in each layer. The edges 
may be in the radial direction, or transversal to it, as illustrated in Figure [M] For 
the network to be critical, and thus uniquely determined by F(7) = F(cr), we must 
nave n odd and / = (n — l)/2 [171 Proposition 2.3, Corollary 9.4], ^ Theorem 2]. 

For the one-sided partial boundary measurements we use a different network 
topology. While conformal or extremal quasiconformal coordinate transformations 
allow for circular networks to be used in the partial measurements case [^, the 
networks with pyramidal graphs F = are more natural [13]. They are shown in 



Figure 15 (left). The pyramidal networks are critical and thus uniquely recoverable 
[T3] . They are natural to use with one-sided partial measurements because the sides 
of the pyramid, where the boundary nodes lie, can be mapped to the accessible 
segment Ba of the boundary. The base of the pyramid consists of interior nodes. 
They model the lack of penetration of the currents in the part of the domain near 
S/, the inaccessible boundary. 

The inversion in the two-sided case is based on two-sided resistor networks [10] 
with graph denoted by F = T^, and n = 2m boundary nodes. There are m nodes on 
each segment of the accessible boundary separated by the leftmost and rightmost 



interior nodes, as illustrated in the right plot in Figure 15 These interior nodes 
model the lack of penetration of the currents in the parts of the domain close to the 
inaccessible boundary. The two-sided network is critical and thus can be recovered 
with e.g. layer peeling [TO] . 
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Figure 15. Resistor networks used for partial boundary mea- 
surements. Left: pyramidal network . Right: two-sided network 
Tn. The boundary nodes Vj, j = 1, . . . , n, n = 10 are numbered 
in circular order and are indicated by x. The interior nodes are 
indicated with o. 



Appendix C. The norm used in the noise model 

Since the boundary B is the unit circle, we associate it with the angle interval 
[0, 27r]. Consider the Fourier series operator Z// : £^ ^ L^[0, 27r], defined by 
(25) 

1 1 /'27r 

{Ufm = ^ V /(fc)e''=^ and its adjoint {U* f) (k) = -= de f{e)e-'''' . 
V2tt V2tt Jo 

The fractional Sobolev norm of / can be written as a weighted £^ norm || ||2, 

(26) \\f\\H^ = \\W-'U*f\\2, where (W'v),^ = {l + k^Y'^ v{k), k G Z. 

The operator norm of a linear operator A : H^^'^{B) H~^^'^{B) is 
(27) 

„,„ WW-^/'Wfh \\W-^/'WAUW-^/'g\\2 

In particular, when A = Ai we have Ai = UKW ^ where Kv{k) = \k\v{k). Thus 

II Al 11^1/2^^:^^-1/2 = 1. 

We approximate the operator norm (27) by the norm || ||, as follows. In an abuse 
of notation, let A^r and Ai be the restrictions of the continuum DtN maps to the N 
uniformly distributed fine grid points on B. Consider the spectral decomposition 

2tt 

(28) --Ai = UT.U\ = 1. 

The approximate norm is given by 

(20) l|AJ= ™p | IK/ + ^')-W(/ + g)-/'.b 

which is equivalent to finding the largest eigenvalue in magnitude of the matrix 
appearing in the numerator above. By construction, we have ||Ai|| = 1. 
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